/*
 * auction_tree.h
 *
 *  Created on: Mar 17, 2011
 *      Author: tqlong
 */

#ifndef AUCTION_TREE_H_
#define AUCTION_TREE_H_

#include <vector>
#include <queue>
#include <algorithm>
#include <limits>

template <typename DATASET>
class TreeDistanceWeight
{
public:
  typedef DATASET                              dataset_type;
  typedef typename DATASET::vec_type           vec_type;
  typedef std::vector<double>                  dvec_type;
protected:
  dataset_type &d1_, &d2_;
  dvec_type *price_;
  double count_;
  double DOUBLE_INFINITY, bigM_;
public:
  TreeDistanceWeight(dataset_type &d1, dataset_type &d2)
  : d1_(d1), d2_(d2), price_(NULL), count_(0),
    DOUBLE_INFINITY(std::numeric_limits<double>::infinity()), bigM_(1e16)
  {
  }

  /** return the weight(i, j) which is the distance of d1[i] and d2[j] */
  double get(int i, int j) {
    count_ += 1;
    return arma::norm(d1_[i]-d2_[j], 2);
  }

  void setPrice(dvec_type& price) { price_ = &price; }
  void setPrice(int i, double p)
  {
    price_->at(i) = p;
  }

  /** return number of points in each set */
  int n1() const { return d1_.n(); }
  int n2() const { return d2_.n(); }
  int n() const { return std::max(n1(), n2()); }

  /** reset the number of accesses to zero */
  void resetCount() { count_ = 0; }
  double count() const { return count_; }

  double cost(int i, int j)
  {
    if (i >= n1() || j >= n2()) return bigM_;
    else return get(i, j);
  }

  void findBest(int i, int &j, int &k, double &v, double &w)
  {
    v = DOUBLE_INFINITY;
    w = DOUBLE_INFINITY;
    for (int x = 0; x < n(); x++) {
      double d = cost(i, x) + price_->at(x);
      if (d < v) {           // if d is smaller than current best
        w = v;               // current best becomes second best
        k = j;
        v = d;               // d becomes current best
        j = x;
      }
      else if (d < w) {      // else, if d is smaller than second best
        w = d;               // d becomes second best
        k = x;
      }
    }
  }
};

template <typename TREECOSTMATRIX>
class AuctionTree
{
public:
  typedef TREECOSTMATRIX                       costmatrix_type;
  typedef std::vector<double>                  dvec_type;
  typedef std::vector<int>                     ivec_type;
protected:
  costmatrix_type &C_;
  ivec_type next_, prev_;
  dvec_type price_;
  std::queue<int> unassign_;
  double DOUBLE_INFINITY, TIGHT_TOLERANCE;
  double bigM_;
  int n_match_;
  double epsilon_;
public:
  AuctionTree(costmatrix_type& C)
  : C_(C), next_(n()), prev_(n()), price_(n()),
    DOUBLE_INFINITY(std::numeric_limits<double>::infinity()),
    TIGHT_TOLERANCE(1e-10)
  {
  }

  int next(int i) const { assert(i < n()); return next_[i]; }
  int prev(int j) const { assert(j < n()); return prev_[j]; }
  int n() const { return std::max(C_.n1(), C_.n2()); }

  void optimize()
  {
    C_.setPrice(price_);
    n_match_ = 0;
    bigM_ = 1e16;
    epsilon_ = 1.0 / n();
    std::fill(next_.begin(), next_.end(), -1);
    std::fill(prev_.begin(), prev_.end(), -1);
    std::fill(price_.begin(), price_.end(), 0.0);
    unassign_ = std::queue<int>();
    for (int i = 0; i < n(); i++) unassign_.push(i);  // maintain an unassign queue
    while (!unassign_.empty()) {
      int i = -1, j = -1, k = -1;
      double v, w;
      i = unassign_.front(); unassign_.pop();         // take 1 unassign item
      findBest(i, j, k, v, w);                        // find best and second best cost(i,j)+price(j)
      C_.setPrice(j, price_[j] + w - v + epsilon_);   // bid new price
      match(i, j);                                    // and match them
    }
  }
protected:
  double cost(int i, int j)
  {
    if (i >= C_.n1() || j >= C_.n2()) return bigM_;
    else return C_.get(i, j);
  }

  void findBest(int i, int &j, int &k, double &v, double &w)
  {
    C_.findBest(i, j, k, v, w);
  }

  void match(int i, int j)
  {
    if (prev_[j] == -1) {
      n_match_++;  // if j is unassign, a new match
      std::cerr << "new match found " << n_match_ << "\n";
    }
    else {
      unassign_.push(prev_[j]);
      next_[prev_[j]] = -1;          // else unassign its current partner
    }
    next_[i] = j;
    prev_[j] = i;
  }
};

#endif /* AUCTION_TREE_H_ */
